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ABSTRACT 

p ^ ■ Low-mass objects embedded in isothermal protoplanetary discs are known to suffer rapid 

inward Type I migration. In non-isothermal discs, recent work has shown that a decreasing 
radial profile of the disc entropy can lead to a strong positive corotation torque which can 
i-C ■ significantly slow down or reverse Type I migration in laminar viscous disc models, depend- 

ing on the amount of viscous and thermal diffusion operating in the planet's horseshoe region. 
Since the latter is a fraction of the pressure scale height of the disc, it is not clear however how 
this picture changes in turbulent disc models. The aim of this study is to examine the impact 
of turbulence on the torque experienced by a low-mass planet embedded in a non-isothermal 
protoplanetary disc. We particularly focus on the role of turbulence on the corotation torque 
whose amplitude depends on the efficiency of diffusion processes in the planet's horseshoe re- 
gion. The main issues we want to address are whether the part of the corotation torque scaling 
with the entropy gradient can remain unsaturated in the presence of turbulence and whether 
the saturation process in non-isothermal discs can be satisfactorily modelled using laminar 
disc models. We performed 2D numerical simulations using a grid-based hydrodynamical 
code in which turbulence is modelled as stochastic forcing. In order to provide estimations for 
the viscous and thermal diffusion coefficients as a function of the amplitude of turbulence, we 
first set up non-isothermal disc models for different values of the amplitude of the turbulent 
forcing. We then include a low-mass planet and determine the evolution of its running time- 
averaged torque. We show that in non-isothermal discs, the entropy-related corotation torque 
can indeed remain unsaturated in the presence of turbulence. For turbulence amplitudes that 
do not strongly affect the disc temperature profile, we find that the running time-averaged 
torque experienced by an embedded protoplanet is in fairly good agreement with laminar disc 
models with appropriate values for the thermal and viscous diffusion coefficients and with the 
formulae of Paardekooper et al. (2011) for the total torque in non-isothermal discs. In discs 
with turbulence driven by stochastic forcing, the corotation torque therefore behaves similarly 
as in laminar viscous discs and can be responsible for significantly slowing down or reversing 
Type I migration. 

Key words: accretion, accretion discs - planets and satellites: formation - hydrodynamics - 
methods: numerical 



1 INTRODUCTION 

In the course of its evolution, the tidal interaction of a planet with 
the protoplanetary disc in which it formed is known to cause a sig- 
nificant change of the planet orbital elements. In particular, low- 
mass planets can experience dramatic changes of their distance 
from the central star under the process of Type I migration (Ward 
1997; Tanaka et al. 2002). Type II migration concerns more mas- 
sive planets that are able to strongly modify the underlying disc 
structure by creating a gap around their orbit (Lin & Papaloizou 
1986). In that case, a planet undergoing Type II migration is locked 
to the disc viscous evolution and migrates inward on a timescale 



corresponding to the disc viscous timescale, provided that the char- 
acteristic disc mass is still larger than the planet mass (Nelson et al. 
2000). 

The gravitational torque exerted by the disc on the planet and caus- 
ing Type I migration consists of two components. The differential 
Lindblad torque results from the angular momentum exchange be- 
tween the planet and the spiral density waves it generates inside the 
disc. For sufficiently low-mass planets, the Lindblad torque is well 
described by a linear analysis, which shows that the specific value 
of this torque scales linearly with disc mass and planet mass and 
as the inverse square of the disc aspect ratio (Tanaka et al. 2002). 
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Although its sign depends on the density and temperature gradients 
inside the disc, the differential Lindblad torque is generally nega- 
tive for typical disc models and is therefore responsible for inward 
migration. 

The corotation torque is due to the torque exerted by the mate- 
rial located in the coorbital region of the planet. A linear analy- 
sis reveals that the corotation torque consists of a barotropic part 
which scales with the vortensity (i.e. the ratio between the verti- 
cal component of the vorticity and the disc surface density) gra- 
dient (Goldreich & Tremaine 1979) plus an entropy-related part 
which scales with the entropy gradient (Baruteau & Masset 2008; 
Paardekooper & Papaloizou 2008). A negative vortensity (resp. en- 
tropy) gradient gives rise to a positive vortensity (resp. entropy) re- 
lated corotation torque. The consequence is that for midly positive 
surface density gradients or negative entropy gradients, a positive 
corotation torque can eventually counteract the effect of a nega- 
tive differential Lindblad torque, which may stall or even reverse 
migration (Masset et al. 2006; Paardekooper & Papaloizou 2009). 
In isothermal discs, the corotation torque is a non-linear process 
generally referred as the horsehoe drag and whose amplitude is 
controlled by advection-diffusion of vortensity inside the horsehoe 
region. In non-isothermal discs, singular production of vortensity 
also arises due to the entropy discontinuity on downstream separa- 
trices (Masset & Casoli 2009; Paardekooper et al. 2010) and in that 
case, the amplitude of the corotation torque is therefore powered 
by advection-diffusion-production of vortensity inside the horse- 
shoe region. 

In the absence of any diffusion processes inside the disc, vortensity 
and entropy gradients across the horseshoe region tend to flatten 
through phase mixing, which causes the two components of the 
horseshoe drag to saturate. Consequently, desaturating the horse- 
shoe drag requires that some amount of viscous and thermal diffu- 
sions are operating inside the horseshoe region. In that case, the am- 
plitude of the horseshoe drag depends on the ratio between the dif- 
fusion timescales and the horseshoe libration timescale and its op- 
timal value, also referred as the fully unsaturated horseshoe drag, is 
obtained when the diffusion timescales are approximately equal to 
half the horseshoe libration time (see e.g. Baruteau & Masset 2012 
for a recent review). In the limit where the diffusion timescales be- 
come shorter than the U-turn timescale, the corotation torque de- 
creases and approaches the value predicted by linear theory. There- 
fore, the corotation torque can be considered as a linear combina- 
tion of the fully unsaturated horseshoe drag and the linear corota- 
tion torque with coefficients depending on the ratio between the dif- 
fusion timescales and the horseshoe libration timescale. Corotation 
torque formulae as a function of viscosity and thermal diffusivity 
were recently proposed by Paardekooper et al. (2011) and Masset 
& Casoli (2010). 

So far, most of the studies aimed at studying the process of satura- 
tion of the corotation torque focused on laminar disc models with 
prescribed values of kinematic viscosity and thermal diffusivity. 
However, it is likely that both angular momentum and heat trans- 
port in protoplanetary discs occur because of turbulence probably 
due to the magneto-rotational instability (MRI, Balbus & Haw- 
ley 1991) and it is not clear how the disc turbulence really im- 
pacts the corotation torque. In the case of isothermal discs, Uribe 
et al. (2011) performed 3D MHD simulations of planet migration 
in stratified discs that suggested the torques coming from the coro- 
tation region can remain unsaturated due to viscous stresses in the 
disc. Baruteau et al. (2011) recently confirmed the existence of an 
unsaturated corotation torque and demonstrated that a horseshoe 
dynamics is indeed at work in turbulent discs with a weak toroidal 



magnetic field. Using 2D hydrodynamical simulations in which 
disc turbulence is modelled through stochastic forcing, Baruteau 
& Lin (2010) showed that in isothermal discs, the running time- 
averaged torque experienced by a protoplanet is in good agreement 
with that obtained from a laminar disc model with a similar vorten- 
sity diffusion coefficient, and in which the total torque consists 
of the Lindblad torque plus the barotropic part of the corotation 
torque. 

In this paper we adopt a similar strategy as in Baruteau & Lin 
(2010) but focus on non-isothermal discs in which the horseshoe 
drag has an additional contribution due to the presence of an en- 
tropy gradient across the horseshoe region. The aim is to investigate 
whether the entropy-related horseshoe drag remains unsaturated 
in presence of turbulence and whether the running time-averaged 
torques in turbulent runs are in good agrement with a laminar vis- 
cous disc model with appropriate values for the vortensity and en- 
tropy diffusion coefficients. 

In the context of planetary population synthesis models aimed at re- 
producing the mass-semi-major axis diagram of exoplanets, it is of 
crucial importance to investigate whether in turbulent disc models 
the entropy-related corotation torque behaves as in laminar viscous 
discs since it may be responsible for significantly slowing down or 
reversing Type I migration in the latter case. 
The paper is organized as follows. In Sect. 2, we give the basic 
equations and notations we use in this work and describe the nu- 
merical set-up in Sect. 3. In Sect. 4, we present the properties of tur- 
bulence in our simulations. In particular, we evaluate the vortensity 
and entropy diffusion coefficients as well as the turbulent Prandtl 
number. In Sect. 5, we estimate the effects of turbulence on the 
entropy-related corotation torque and compare with the torque for- 
mulae of Paardekooper et al. (2011) and with laminar calculations 
as well. We finally discuss our results and draw our conclusions in 
Sect.6. 



2 PHYSICAL MODEL 
2.1 Basic equations 

We adopt a 2D adiabatic disc model for which all the physical quan- 
tities are vertically averaged. We work in a frame rotating with an- 
gular velocity Q, p , and adopt cylindrical polar coordinates (R,(/>) 
with the origin located at the position of the central star. The conti- 
nuity, equation of motion and energy equation respectively read: 

g+V-(Z?) = 0, (1) 
at 

dv Vp 

_+(v?.V)i?=-^-VO, (2) 

and, 

dS 

— + (^V)S=0 (3) 
at 

where 2 is the disc surface density, v the velocity, P the pressure 
and S the entropy which is defined as: 



where y ad is the adiabatic index. Here we use an ideal gas equation 
of state p = ZR g T/fi where R g is the gas constant, T the temperature 
and ji the mean molecular weight. 
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In Eq.[2] O is the gravitational potential and reads: 



<t> : 



GM+ 1 , , 



(5) 



where M t is the mass of the central star and CD;,,,/ is an indirect 
term which accounts from the fact that the frame centered on the 
star is non-inertial (e.g. Nelson et al. 2000). Effects resulting from 
turbulence are modelled by applying the external turbulent poten- 
tial (Dftjrf, of Laughlin et al. (2004) (see Sect. 3). Previous studies 
(e.g. Baruteau & Lin 2010; Pierens et al. 2011) have shown that 
the perturbations induced by this stochastic forcing term have sim- 
ilar statistical properties to those resulting from 3D MHD simula- 
tions. In isothermal discs, it is well known that fluctuations arising 
from turbulence generate transport of angular momentum through 
Reynolds and Maxwell stresses (e.g. Papaloizou & Nelson 2003). 
For non-isothermal discs, additional turbulent energy transport is 
expected as a consequence of entropy perturbations. 



2.2 Averaged quantities and turbulent diffusion coefficients 

To estimate the rate of angular momentum and energy transport 
resulting from turbulence, we adopt a common procedure and make 
use of averaged quantities. The mass-weighted azimuthal average 
of the quantity Q(r, <p, t) is given by: 



Q(R,t) 



J E(fl, (p,t)Q(R,if>,t)d<f> 
f L(R, (p, t)d(p 



(6) 



where the integral is performed over 2n in azimuth. Setting vr = 
vr + Svr and = + Sv^, where 6v R and cSv^ are the radial and 
azimuthal velocity fluctuations respectively, and averaging the az- 
imuthal component of the equation of motion (Eq. |2) gives the fol- 
lowing equation for the conservation of the specific angular mo- 
mentum j = Rv$ (e.g. Fromang & Nelson 2006): 



d 



1 d 



1 d 



_( S J )+ __ ( ^ J) = ___(^, V>) ( 7) 

The previous equation can be re-written as (Balbus & Papaloizou 
1999): 

1 d 



RdR 



(8) 



where a R is the Reynolds stress parameter which is defined as: 



a R 



(9) 



Adopting a similar averaging procedure for the energy equa- 
tion leads to the following relation for the conservation of entropy: 



arising from the turbulence (see Sect.[3}. Turbulent heating remains 
small, however, for the turbulence amplitudes considered in this 
study. To estimate the effective thermal diffusion coefficient in the 
simulations, we therefore simply neglect the turbulent heating term 
and make use of Eg. 1121 



3 NUMERICAL SET-UP 

Simulations were performed with the FARGO and GENESIS nu- 
merical codes which solve the equations for the disc on a polar 
grid. These codes employ an advection scheme based on the mono- 
tonic transport algorithm (Van Leer 1977) and include the FARGO 
algorithm (Masset 2000) to avoid time step limitation due to the 
Keplerian orbital velocity at the inner edge of the grid. The energy 
equation that is implemented in both codes reads: 



— + V • (ev) 
dt 



-(.7a 



(13) 



where e is the thermal energy density and Ql ulk is the heating 
source term by artificial viscosity and provided by shocks arising 
from turbulent fluctuations. This term is handled by adopting a 
standard von Neumann-Richtmyer artificial bulk viscosity (see 
Stone & Norman 1992) with Ci = 1-4, where Co represents the 
number of zones over which artificial viscosity will spread a shock. 
Here, no cooling term is included as an additional source term, 
which means that reaching a thermal equilibrium for the disc is not 
possible and that the temperature profile may be progressively al- 
tered by turbulent heating, especially in the case of strong turbulent 
fluctuations. However, from the results of simulations presented in 
Sect. [5] we anticipate that using a cooling term is not necessary 
since for most of the turbulence amplitudes considered in this 
study, turbulence heating has little impact on the temperature 
profile over the runtime covered by the calculations. In fact, not 
considering such a term enables to investigate effects resulting 
from turbulent diffusion only without continuously changing the 
temperature profile close to the planet. We further note that the 
results presented in Sect.[5]may depend strongly on the expression 
employed for the cooling term. 

The computational units that we adopt are such that the unit of 
mass is the central mass M,, the unit of distance is the initial orbital 
radius R p of the planet and the unit of time is fl" 1 = (GM,/R 3 p )~ 1/2 . 
In the simulations presented here, we use Nr = 512 radial grid 
cells uniformly distributed between Ri„ = 0.4 and R m „ = 1.8, and 
N$ = 1200 azimuthal grid cells. Wave-killing zones are employed 
for R < 0.5 and R > 1.58 to avoid wave reflections at the disc 
edges (de Val-Borro et al. 2006). 



|(^) + i|(^5) = -I|( W 5) (10) 

which can be recast as: 

d - 1 d -Id dS 

Jt^^RdR^^RM^^ 

where k is the diffusion coefficient associated with the entropy and 

is given by: 



6v R 6S 
k = — 



dS/dR 



(12) 



Here, we note that the previous equation was derived assuming 
an adiabatic disc model. This neglects turbulent heating, which is 
modelled in our simulations by artificial viscous heating at shocks 



The initial disc surface density profile is S = S p x {R/Rp)^ 
with = 5 x 10~ 4 and cr = 1.5 by default. For such a disc surface 
density profile, the initial vortensity-related part of the corotation 
torque cancels out so that the corotation torque consists only of its 
entropy-related part. The initial disc aspect ratio is h = h p x (R/R P Y 
with h p = 0.03 and where / is the disc flaring index which is set to 
/ = 0.3 in this work. This corresponds to an initial temperature pro- 
file decreasing as R~P with ft = 1 - 2/ = 0.4. For an adiabatic index 
7 = 5/3 which corresponds to the value adopted here, the initial en- 
tropy profile therefore varies as R ^ with £ = [}-(y- \)cr = -0.59. 
As will be justified in Sect.|5] such a positive gradient for the initial 
entropy profile is chosen to minimize the maximum convergence 
time of the turbulent runs, since both the entropy-related corotation 
torque and the differential Lindblad torque are negative in that case. 
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In the simulations presented here, turbulence is modelled by 
applying at each time-step a turbulent potential O niri , to the disc 
(Laughlin et al. 2004, Baruteau & Lin 2010) and corresponding to 
the superposition of 50 wave-like modes. This reads as: 

50 

<t>, urb (R, fa t) = jR 2 Q 2 J] A k (R, fa t), (14) 

k=l 

with 

A* = &e "* cos(m k (j> - fa - Q. k t k ) sm(nf k /At k ). (15) 

In Eq.Q/J] £ k is a dimensionless constant parameter randomly 
sampled with a Gaussian distribution of unit width. R k and fa 
are, respectively, the radial and azimuthal initial coordinates of the 
mode with wavenumber m k , <r k = nR k /4m k is the radial extent of 
that mode, and Q. k denotes the Keplerian angular velocity at R = R k . 
Both R k and fa are randomly sampled with a uniform distribution, 
whereas m k is randomly sampled with a logarithmic distribution 
between m k = 1 and m k = 96. Each mode of wavenumber m k starts 
at time t = to tk and terminates when f k = 1 - > At k , where 
At k = 0.2nR k /m k c s denotes the lifetime of mode with wavenumber 
m k , c s being the sound speed. Such a value for Atk yields an auto- 
correlation timescale r c ~ 0.5T orb , where T orb is the orbital period 
at R = 1 (Baruteau & Lin 2010, and see Sect.gl}. 
Following Ogihara et al. (2007), we set At = if m k > 6 to save 
computing time. As noticed by Baruteau & Lin (2010), such an 
assumption is supported by the fact that a turbulent mode with 
wavenumber m has an amplitude decreasing as exp(-m 2 ) and a life- 
time oc 1/m, so that the contribution to the turbulent potential of a 
high wavenumber turbulent mode is relatively weak. 
In Eq. [14] y denotes the value of the turbulent forcing parame- 
ter, which controls the amplitude of the stochastic density pertur- 
bations. In the simulations presented here, we used values for y 
ranging from 10~ 5 to 2 x 10~ 4 . Given that y is related to the dimen- 
sionless a viscosity of Shakura & Sunyaev (1973) by the relation 
a ~ I20(y/h p ) 2 (see Sect. |4.1K the latter values for y correspond to 
a ranging from ~ 10~ 5 to ~ 5 x 10~ 3 . 

In the following, we first estimate how both the viscous and ther- 
mal diffusion coefficients induced by turbulence depend on the am- 
plitude of the turbulent potential y. We then present the results of 
calculations of a low-mass planet embedded in a non-isothermal 
disc subject to stochastic forcing. We consider different values for 
y and study how the running time-averaged torque deduced from 
our simulations compare with laminar runs and with the analytical 
formulae of Paardekooper et al. (201 1) for the torque due to Lind- 
blad resonances and horseshoe drag in the presence of viscous and 
thermal diffusion. 



4 TURBULENCE PROPERTIES 

4.1 Estimation of the turbulent thermal diffusivity 

In non-isothermal discs, the main effect of the applied turbulent 
potential is to generate velocity and temperature fluctuations 
responsible for both entropy and angular momentum transport. 
In a previous study dedicated to isothermal discs, Baruteau & 
Lin (2010) showed that for the adopted turbulent potential, the 
resulting Reynolds stress parameter scales as y 2 . Here, in order to 
investigate the dependence of the entropy diffusion coefficient on 
the amplitude of turbulence y, we performed a series of simulations 
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Figure 1. Upper panel: Reynolds stress parameter a(R) as a function of 
radius for y = 8 X 10~ 5 and time-averaged over 700 T orb . Lower panel: 
same but for the entropy diffusion coefficient k{R) defined by Eq . 1121 

in the non-rotating frame of a non-isothermal disc with turbu- 
lence driven by stochastic forcing, varying the value for y from 
7 = 2 x 10~ 5 to y = 2 x 10~ 4 and using the disc model described in 
Sect. [3] From the results of these simulations, the Reynolds stress 
parameter a R (R,t) and the entropy diffusion coefficient K(R,f) are 
computed using Eqs.|9]and[T2]respectively. To examine how these 
two quantities vary with radius, a time average over 700 T orb is 
performed. For y = 8 x 10~ 5 , Fig. \T\ displays the time-averaged 
coefficients a R (R) and k(R) as a function of radius. Clearly, both 
quantities appear to be approximately uniform with radius, except 
in the vicinity of the wave-killing zones. 

We now consider an additional space average of a R (R) and 
k(R) in the range 0.9 < R < 1.1 and in the following, we simply 
denote by a R and k the resulting numbers. In simulations with an 
embedded planet presented in Sect.|5] the planet is on a fixed circu- 
lar orbit located at R p = 1. Given that the corotation torque expe- 
rienced by the planet depends on the viscous and thermal diffusion 
timescales across the horseshoe region whose half-width is a frac- 
tion of the disc scale height H for the planet masses we consider, 
the range over which the space average is performed therefore cor- 
responds to the disc region which contributes the most to the torque 
exerted on the planet. The upper panel of Fig.|2]represents a R as a 
function of y. As expected, a R appears to scale as y 2 and is found 
to be well approximated by: 

<*r *Mt~\ > ( 16 > 
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Figure 2. Upper panel: Averaged Reynolds stress parameter a as a func- 
tion of turbulence amplitude y. Lower panel: Averaged entropy diffusion 
coefficient a: as a function of turbulence amplitude y. 



which is in excellent agreement with the relation derived by 
Baruteau & Lin (2010). A slight discrepancy arises because we 
average with radius in the range 0.9 < R < 1.1 whereas the 
space average is performed over the disc region located outside 
the wave-killing zones in Baruteau & Lin (2010). While in laminar 
discs vortensity diffusion is controlled by the kinematic viscosity, 
this is not necessary the case in turbulent discs where the vorten- 
sity's diffusion coefficient D can differ from the turbulent viscos- 
ity v R = a R c s H associated with the turbulent diffusion of angular 
momentum. For instance, Baruteau & Lin (2010) have shown that 
in isothermal discs, the ratio S c = v R /D usually referred as the 
Schmidt number is S c ~ 0.25. Although not shown here, we have 
checked that this value also holds in non-isothermal discs and in 
the following, we will make use of the dimensionless vortensity's 
diffusion coefficient that we define as a D = S^ l a R ~ 4a R . Given 
that in MHD calculations, the total viscous stress parameter corre- 
sponding to the dimensionless a viscosity of Shakura & Sunyaev 
(1973) is ~ 4a R (Fromang & Nelson 2009), this suggests that our 
value for a D is close to that corresponding to the standard a vis- 
cous stress parameter. Using Eq.QJo] the dimensionless viscosity of 
Shakura & Sunyaev (1973) can consequently be well approximated 
by: 



(i 1201 — 



(17) 



We plot the entropy diffusion coefficient k as a function of y in 
the lower panel of Fig. [2] Again, k is found to increase as y 2 , which 
is in agreement with the expectation that both 6v r and SS scale with 
y. The best-fitting solution is found to be: 
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Figure 3. Perturbed entropy profile at / = 4.5T or i, averaged over six differ- 
ent realizations with y = 6 X 10~ s (left panel) and y = 10~ 4 (right panel) 
and in which a 10% negative perturbation is applied to the background en- 
tropy profile at Rj = I. The solid line corresponds to the analytical fit given 
by Eq. |19| with k given by Eq. |18| 
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Figure 4. Upper panel: Surface density profile at t = 700 T or j, for the 
different values of y we considered. Lower panel: Same but for the entropy 
profile. 



To check the consistency of the previous estimation, we fol- 
lowed the approach of Baruteau & Lin (2010) and performed a 
series of runs varying the value for y and in which a 10% nega- 
tive perturbation with a (S-function profile is applied to the initial 
entropy profile at Rj ~ 1. This is achieved by applying a positive 
perturbation to the density profile together with a negative temper- 
ature perturbation with same amplitude, in such a way that the ini- 
tial pressure perturbation is zero. The evolution of the perturbed 
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entropy profile is then compared with the function: 

^ 

5S(R,t) 



6S AR l-(R-R d Y 
, expf 



4kt 



(19) 



where 5S o is the initial entropy perturbation at R = R d , AR is the 
radial resolution of the grid and k is the estimation of the entropy 
diffusion coefficient given by Eq. [T8] The results of this set of sim- 
ulations are presented in Fig. [3] where we plot for y = 6 x 1(T 5 
and y = 1CT 4 both the perturbed entropy profile averaged over six 
different realizations at t = 4.5 T or t and the fitting formula given by 
Eq. [19] Clearly, the perturbed entropy profile closely matches the 
expected one, which indicates that Eq. [T8] provides a fairly good 
estimation of the effective entropy diffusion coefficient. 
It is worthwhile to notice that this model can suffer from the time 
evolution of the entropy profile due to artificial viscous heating at 
shocks. This is illustrated in Fig. [4] which shows the surface den- 
sity and entropy profiles at r = 700 T m ], for the different values of 
y we considered. Although the surface density profile remains al- 
most constant, we see that the entropy profile can be significantly 
altered in turbulent runs with high values of y. This is due not only 
to turbulent diffusion which makes the entropy profile steeper but 
also to turbulent heating which tends to increase the local value for 
the entropy. Looking at Fig. [4] we see that effects resulting from 
turbulent heating are more pronounced near the radial boundaries. 
However, we are confident that this process does not strongly affect 
our estimations of a and k since these are determined using space 
averages in the range 0.9 < R < 1.1, which corresponds to disc 
regions where the entropy profile is only slightly modified over the 
runtime of the simulations. To check this, we computed these two 
quantities using normalizations by the initial pressure or entropy 
profiles and found good agreement with the estimations given by 
Eqs.ll7landll8lfor which normalizations by time averaged profiles 
were considered. 



4.2 Determination of the turbulent Prandtl number 

In this section, we investigate how the turbulent Prandtl number P r 
depends on the amplitude of the turbulent forcing y. The turbulent 
Prandtl number is commonly defined as: 



Pr 



Vr 

K 



a R c s H 



(20) 



Since both a R and k scale with y 2 , we expect P r to not depend 
on the value for y. This is confirmed by inspecting Fig. [5] where 
P r is displayed as as a function of y and which reveals that the 
turbulent Prandtl number is P, ~ 0.3 on average. Such a value is 
close to that derived by Ruediger et al. (1988) who estimate that 
a value P r ~ 0.1 is necessary to explain the outburst behaviour 
of cataclysmic variables. In the case of shear-driven hydrodynamic 
turbulence, it is also of interest to mention that Lathrop et al. (1992) 
measured a turbulent Prandtl number P, ~ 0.176 in experiments of 
Couette-Taylor flows at Reynolds numbers well beyond the onset 
of chaos. 



4.3 Turbulent transport of entropy in a radiative 
protoplanetary disc 

As mentioned above, no cooling source term is included in the en- 
ergy equation that is employed in this work. However, in a more 
realistic disc model in which turbulent effects and radiative pro- 
cesses are both included, we expect a thermal equilibrium state to 
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Figure 5. Estimated turbulent Prandtl number as a function of turbulence 
amplitude y. 



be reached when radiative cooling compensates turbulent heating. 
For such a radiative turbulent disc, diffusion of entropy is due not 
only to turbulent transport as it is the case for the idealized disc 
considered here but also to radiative diffusion. For an active lami- 
nar disc model with turbulent heating balanced by radiative losses, 
the thermal diffusion coefficient K ra a associated with radiative trans- 
port can be expressed as (Paardekooper et al. 201 1): 



, 2V3 2 
Kmd = 9y ad (y ad - l)v| 1 + + 



(21) 



where r is the optical depth and v = ac s H the kinematic viscosity. 
Noting that v = S~ l P r k from the results of the previous section 
and using Eq. [2J] we obtain the following expression for the ratio 
between the thermal diffusion coefficient associated with radiative 
effects and that associated with the turbulent transport of entropy: 

K -f = 9S-JP r y ad (y ad - 1) (l + ^ + Jjj (22) 

Previous equation shows than in optically thick inner regions of 
protoplanetary discs, turbulent and radiative diffusion coefficients 
will tend to be of the same order of amplitude whereas in optically 
thin outer regions, radiative transport will be much more efficient 
than turbulent transport for diffusing entropy. 



5 EFFECTS OF TURBULENCE ON THE 

ENTROPY-RELATED COROTATION TORQUE 

In this section, we evaluate the effect of turbulence on the entropy- 
related corotation torque. To investigate this issue, we performed 
simulations of a low-mass planet on a circular orbit at R p = 1 and 
(p p = n and embedded in a non-isothermal disc subject to stochastic 
forcing. We used values for the amplitude of the turbulent potential 
ranging from y = 10~ 5 to y = 2 x 10~ 4 . The disc model is the 
same as described in Sect. [3] with power-law indexes for the ini- 
tial surface density, temperature and entropy profiles of a = 1.5, 
P = 0.4 and £ = -0.59 respectively and an aspect ratio at the 
planet position of h„ = 0.03. We considered planets with mass ra- 
tio q = 5 x 10~ 6 and q = 10~ 5 which corresponds to planets with 
masses of 1.6 M e and 3.2 M e if the central star has a solar mass. 
For q = 5 x 10~ 6 and for a softening parameter b = 0.4h p which is 
the value set in this work, the half-width of the horseshoe region is 



0.013 (Paardekooper et al. 2010) which 



means that the horseshoe region is resolved by about 12 grid cells 
in the radial direction. Also, we note that in calculating the torque 
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experienced by the planet, we include the disc material located in- 
side the Hill sphere of the planet. 

5.1 Total torque exerted on the planet 

As shown by Baruteau & Lin (2010) in the case of isothermal discs, 
turbulence can unsaturate the vortensity-related part of the corota- 
tion torque in the same way as an alpha viscosity does in a laminar 
disc. In the following, we make the assumption that thermal dif- 
fusion resulting from turbulence can also unsaturate the entropy- 
related part of the corotation torque and assume that this effect 
can be modelled by an entropy diffusion coefficient. This will be 
checked in the next section. 

In the presence of turbulence, the torque exerted on a low-mass 
planet embedded in a turbulent disc can be considered as the sum of 
the differential Lindblad torque T L , the corotation torque T c and the 
turbulent torque Y turB . In the adiabatic limit and for a smoothing pa- 
rameter b = 0.4h p , the Lindblad torque is given by (Paardekooper 
et al. 2010): 



where To is given by: 



(-2.5- 1.70 + O.lcr) — 

Jad 



(23) 



(24) 



In the previous equation, the subscript "p" indicates that quantities 
are evaluated at the orbital position of the planet R p . 
Provided the surface density profile is not strongly affected by tur- 
bulence, we expect the barotropic part of the corotation torque to 
cancel for our disc model in which cr = 1.5, since the latter scales 
with the vortensity gradient. In that case, the corotation torque con- 
sists only of its entropy-related part T c e and the total torque T ex- 
erted on the planet simplifies into: 



r = r L + r„ + r„ 



(25) 



For diffusion timescales across the horseshoe region shorter than 
the horseshoe libration period Tm, = &nR p /3£l p x s but longer than 
the horseshoe U-turn timescale rp_,„ r „ ~ h p Tm, (Baruteau & Mas- 
set 2008), the entropy-related horseshoe drag is close to its fully 
unsaturated value given (Paardekooper et al. 2010): 



T/j.y,e 



(26) 



r tm\ 

Jad \ Jad I 

while for diffusion timescales smaller than the U-turn timescale, 
T c e decreases toward the entropy-related part of the linear corota- 



tion torque T c j ilLe with: 



r c , iin , e = — 2.2-— u 

Jad \ 7ad) 



(27) 



A model of the entropy-related corotation torque including its de- 
pendence on the diffusion timescale can be obtained by writing 
T ce as a linear combination of and T c j inf . Paardekooper et al. 
(201 1) have shown that T ce can be indeed well approximated by: 



(28) 



r c , f = F(p v )F(p x ) ^G(p Y )G{p x )T,, SJ , 
+ ^/(l - K( Pv ))(l - K(p x )T cJw , e 

where p v is the saturation parameter related to viscosity and is given 
by: 



27TV 



(29) 



and p x is the saturation parameter related to thermal diffusion with: 



2 mprf 



(30) 



where v = ac s H with a given by Eq. [TTJand where k is related to 
the amplitude of turbulence by Eq.[T8] In Eq. 28, the function F(p) 
of the parameter p is defined by: 

1 



p 1 + GV1.3) 2 

while the functions G{p) and K(p) are given by: 
G(p) 



4/3 



if p < V8/(45?r) 



,-8/3 



if p > V8/(45tt) 



and: 



J 25 V 28 / J 



} 3/2 , 



if p < V28/(45tt) 



ifp> V28/(45w) 



(31) 



(32) 



(33) 



Regarding the turbulent torque r,„,;,; its amplitude is given by (Nel- 
son 2005; Baruteau & Lin 2010): 



IT^uri?! — 0"turb 



-1/2 



(34) 



where cr nirb ~ 2.4 x 10 2 S p #yfi£n 2 (Baruteau & Lin 2010) 
is the mean deviation of the turbulent torque distribution and t c ~ 
0.5 T or i, is the correlation time of the turbulence. Using the previous 
equation leads to an estimation of the convergence time t conv of the 
simulations over which \r turc ,\ is a small fraction / of \T L + r c>e |, 
Using Eq.(34] the convergence time is given by: 

v 2 



& turb 

T L + fZ 



(35) 



Previous equation shows that the maximum convergent time of the 
turbulent runs is reduced in the case where both the entropy-related 
corotation torque and the differential Lindblad torque have the same 
sign. In order to minimize t com , and given that the differential Lind- 
blad torque is negative, we employed a positive initial entropy gra- 
dient so that the entropy-related corotation torque is negative as 
well. Using the expression for T L given by Eq. [23] and taking the 
fully unsaturated value for the entropy-related corotation torque 
given by Eq. [26] we estimate the convergence time to be approxi- 
mately given by: 

,2 



0.2 - T, 



(36) 



where we assumed / = 0.1. For q - 5 x 10 -6 and y = 10~ 4 , this 
gives t com ~ 100 T or t Over timescales typical to those covered 
by our simulations (~ 1500 T orb ) and provided that there is no 
significant coupling between the Lindlad torque and the turbulent 
torque, we therefore expect the steady-state torque experienced by 
the planet to be approximately given by: 



r ~ r L + r c>e 



(37) 



where T L and T c e are respectively given by Eqs. [23] and [29] We 
emphasize that the previous equation is valid in the case where 
the effect of turbulent thermal diffusion on the desaturation of the 
entropy-related corotation torque is similar to that corresponding 
to a laminar entropy diffusion coefficient. The aim of the following 
section is to check this assumption. 
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5.2 Time evolution of the running time-averaged torque 

For q = 5 x 10~ 6 , the time evolution of the running time-averaged 
torque over 1500 orbits for turbulent runs with y in the range 
10~ 5 - 2 x 10~ 4 is displayed in the left panel of Fig. [6] From Eq. 
1361 we expect that the running time-averaged torques have con- 
verged over this timescale except for the run with y = 2 x 10" 4 in 
which the temperature profile is found to be continuously altered 
by turbulence (see Sect. [53). Also shown is the result of an invis- 
cid simulation with y = in which the torque saturates toward the 
differential Lindblad torque T L . In that case, T L is found to slowly 
increase by ~ 2% between 1000 and 1500 orbits due to the con- 
tinuous formation of a gap around the planet. We note that such a 
process was also observed in the isothermal turbulent simulations 
ofBaruteau&Lin (2010). 

In turbulent runs, the running time-averaged torque typically 
reaches a quasi-stationary value after ~ 200 orbits for moderate 
values of y, which is relatively consistent with the estimation of the 
convergence time given by Eq. [36] Up to y = 10~ 4 , we see that 
this stationary value decreases as y increases whereas for higher 
values of y, it tends to increase with y. As the differential Lindblad 
torque should not be significantly altered by turbulence (Baruteau 
& Lin 2010, Baruteau et al. 2011) and since the barotropic part 
of the corotation torque is expected to cancel for the chosen surface 
density profile (cr = 1.5), this seems to suggest that turbulence does 
unsaturate the entropy-related corotation torque. In that case, we in- 
deed expect that as y increases, the torque experienced by the planet 
decreases from the differential Lindblad torque down to value close 
to that corresponding to the sum of the differential Lindblad torque 
plus the unsaturated entropy-related horseshoe drag. Such a value 
is reached when the thermal diffusion timescale across the horse- 
shoe region r d = x 2 /k is approximately equal to half the libration 
timescale T;«, = SnR p /(3n p x s ). Using Eq. [T8] we estimate the un- 
saturated entropy-related horseshoe drag to be obtained for: 

/ \ 3/4 

7-0.07^) (38) 

Using Eq.QT] this corresponds to: 

a ~ 0.6q m h- 712 (39) 

For q = 5 x 10~ 6 , this gives y ~ 1.1 X 10~ 4 (a ~ 1.6 x 10~ 3 ) 
which is close to the value observed in the simulations. For higher 



values of y, the diffusion timescale across the horseshoe region can 
approach the U-turn timescale Tu-tum ~ h p Tiu, (Baruteau & Masset 
2008) and the amplitude of the entropy-related corotation torque 
starts to decrease toward a value corresponding to that predicted by 
the linear theory. 

The results for q = 10~ 5 are shown in the right panel of Fig.[6]which 
presents the time evolution of the running time-averaged torque for 
y < 2 x 10~ 4 . For y = 0, it appears that the torque oscillates around 
a value corresponding to the differential Lindblad torque. As dis- 
cussed in Baruteau & Lin (2010), these oscillations are due to the 
formation of vortices flowing along the edge of the gap created 
by the planet (Li et al. 2005). However, the running time-averaged 
torque does converge toward the differential Lindblad torque but 
compared with the case q = 5 x 10~ 6 , the torque saturates much 
more slowly. For y ± 0, we estimate using Eq.[38]that the value of 
7 from which the entropy-related horseshoe drag starts to decrease 
toward its linear value is y ~ 1.8 x 10~ 4 . This is in decent agree- 
ment with the results of the simulations for q = 10~ 5 which show, 
in the limit of y < 2 x 10~ 4 , a tendency for the torque amplitude to 
increase with y. Again, this enlightens the process of desaturation 
of the entropy-related corotation torque by turbulence. 



5.3 Comparison with laminar viscous disc models 

In this section, we compare the values for the running time- 
averaged torques deduced from the turbulent runs with the 
torque formulae for non-isothermal Type I migration derived by 
Paardekooper et al. (2011) and with laminar simulations with con- 
stant values for the kinematic viscosity and thermal diffusivity. In 
Fig.[7]are plotted, for q = 5 x 10~ 6 (left panel) and q = 10~ 5 (right 
panel), the running time-averaged torques computed at t = 1500 
r or 2, as a function of the estimated entropy turbulent diffusion co- 
efficient k. The top axis gives the corresponding value for the y 
parameter which is related to k by the relationship given in Eq. 
n~8l The solid line in Fig. [7] depicts the approximation for the to- 
tal torque given by Eq. [37] where we remind the reader that T L is 
the Lindblad torque and T c e the entropy-related corotation torque 
given by Eq. [29] In this latter relation, we computed the satura- 
tion parameter related to viscosity p v by setting v = S^P r K ~ 1.2k. 
Moreover, we used the results from the run with 7 = and in which 
the corotation torque is expected to saturate to determine T L . From 
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Figure 7. Running time-averaged specific torque (stars) in turbulent runs as a function of the turbulence amplitude y (top axis) and as a function of the effective 
entropy diffusion coefficient k (bottom axis) computed using Eq. [18] Diamonds correspond to the results of laminar simulations with values for the viscous 
and thermal diffusion coefficients equal to those effective in turbulent runs. The initial surface density and entropy profiles in laminar runs correspond to the 
profiles in turbulent runs time averaged over 1500 T or b- The solid line depicts the analytic estimation of Paardekooper et al. (201 1) expected for the torque as 
a function of the thermal diffusivity k, for a kinematic viscosity v = S c P r K = 1.2k and using the initial surface density and entropy profiles in turbulent runs. 
The upper and lower dashed lines depict the values of the differential Lindblad torque and the fully unsaturated torque respectively. 



this simulation we measure T L ~ 1.54r for q = 5 x 1CT 6 and 
T L ~ 1.34ro for q = 1CT 5 which differs from the estimation for the 
Lindblad torque given in Eq.[23]by ~ 15% and ~ 25% respectively. 
We interpret this slight discrepancy as being due to a change of the 
temperature profile in the vicinity of the planet. This is exempli- 
fied in Fig.[8]which displays for the two values of q we considered 
the temperature profile in the case where y = 0. It reveals that the 
temperature tends to increase at distances ±2H from the planet, cor- 
responding to a local increase of the disc aspect ratio and resulting 
consequently in a decrease of the Lindblad torque. 
For the run with y = 0, the fully unsaturated horseshoe drag was 
evaluated to be r,,.„ ~ 1.82r for q = 5x 1CT 6 and r,,„ = 1.92r 
for q = 1(T 5 which are values close to the estimation given by Eq. 
[26] which predicts ~ 1.7Fo. Regarding the linear corotation 
torque, it can be crudely estimated by subtracting the differential 
Lindblad torque from the total torque at t ~ 2 T orB ( that is, t ~ 
0.5r f y_„„„; Baruteau & Masset 2008; Paardekooper & Papaloizou 
2008; Casoli & Masset 2010) and this gives T cM<e ~ 0.7ir for 
q = 5 x lO" 6 and T cMe ~ O.S6T a for q = 10~ 5 while using Eq.l27l 
leads to T c j jn e ~ 0.46r . 

For modest values of y, we see that good agreement is obtained 
between the results of the turbulent runs and the analytic formulae 
of Paardekooper et al. (201 1). The basic reason for the differences 
at higher values of y possibly arises from the evolution of the disc 
temperature profile. This is illustrated in Fig. [9] which shows for 
q = 5 x 10~ 6 the entropy and surface density profiles time-averaged 
over 1500 T orB . We note that in Fig. [9] the entropy profile corre- 
sponding to y = 2 x 10~ 4 was offset by a factor of 2 for clarity. 
Although the density remains almost unchanged, the entropy pro- 
file can be significantly affected in the vicinity of the planet for 
y > 10~ 4 and this should in principle affect the value for the to- 
tal torque. For y = 10~ 4 , we find that with respect to an inviscid 
run, the entropy gradient is increased by ~ 20% at the location 
of the planet, which should increase the amplitude of the (nega- 
tive) entropy-related corotation torque by the same amount. The 
corresponding change on the slope of the temperature profile is es- 
timated to lead to a ~ 13% decrease in the amplitude of the (nega- 
tive) differential Lindblad torque. Interestingly, such changes in the 



values for the corotation and the differential Lindblad torques are 
found to almost cancel each other which suggests that in that case 
the modification of the temperature by turbulence has little impact 
on the total torque experienced by the planet. In the calculation with 
7 = 2 x 10~ 4 , the surface density at R = 1 in only slightly reduced 
by ~ 3% with respect to the initial profile but the temperature is 
found to be higher by a factor of ~ 2, with the consequence that the 
torque exerted on the planet is weaker by the same factor compared 
to a laminar run. Moreover, we quantify the corotation torque and 
the differential Lindblad torque to be altered by ~ 57% and ~ 65% 
respectively due to the modification of the temperature slope at the 
position of the planet. From these values, we estimate that for the 
run with y = 2 x 10~ 4 and q = 5 x 10~ 6 , the specific torque exerted 
on the planet would be ~ 10~ 5 in the absence of modification of the 
temperature profile by turbulence, which is somewhat larger than 
the value predicted by Eq.[37] We note that differences between the 
time-averaged turbulent torques and the analytical estimations may 
be also related to uncertainties in Eqs.[T7] and \TE\ which relate the 
value of y to those of a and k. 

We now turn to the question of how the torque dependence 
upon turbulence strength in turbulent simulations compares with 
the torque dependence upon kinematic viscosity and thermal diffu- 
sivity in laminar runs. To investigate this issue, we have performed 
an additional suite of laminar simulations in which we vary the 
value for the thermal diffusivity while the kinematic viscosity is set 
to v = S~ [ P r K ~ 1.2k. The values for the thermal diffusivity are 
chosen in such a way as to correspond to the effective values for 
thermal diffusivity in turbulent runs and estimated using Eq.[T8] In 
order to take into account the alteration of the temperature profile at 
high levels of turbulence, the inital surface density and temperature 
profiles in these laminar simulations consist of the surface density 
and temperature profiles in turbulent runs time-averaged over 1500 
orbits. The laminar torques typically reach a quasi-stationary value 
by 500 To,!; and they are directly compared to the running time- 
averaged turbulent torques at t = 1500 T orb in Fig. [7] Again, for 
the two values of q we consider, relatively good agreement is ob- 
tained between the laminar and turbulent runs, even for high values 
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Figure 10. Left panel: Snapshot of the disc perturbed entropy at t = 800 T or b for a laminar run with q = 5 X 10 -6 and with constant thermal and viscous 
diffusion coefficents with values k = 4.5 X 10~ 7 and v = 5.4 X 10~ 7 respectively. Streamlines are overplotted as black lines. Middle panel: same but for a 
turbulent run with y = 6 X 10~ 5 . Right panel: Contour plot of the perturbed entropy for the same turbulent run but time-averaged between 750 and 850 orbits 
and over 6 different realizations. Averaged streamlines are overplotted as black lines. 
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Figure 8. Temperature profile at t = 1000 T orb for q = 5 X 10~ 6 black 
line and q = 10~ 5 (red line). The dashed line corresponds to the initial 
temperature profile. 



of the thermal diffusivity for which the temperature profile can be 
strongly modified in turbulent runs. One can notice however a ten- 
dency for the turbulent torques amplitudes to be slightly smaller 
than the laminar torques, which would suggest that for similar val- 
ues of the thermal and viscous diffusion coefficients, the corotation 
torque is slightly more saturated in turbulents discs than in laminar 
viscous disc models. We note that such a trend was not observed 
in previous studies of the vorticity-related corotation torque in tur- 
bulent discs (Baruteau & Lin 2010; Baruteau et al. 2011) which 
suggested that the turbulent torque experienced by a protoplanet 
embedded in a isothermal turbulent disc is in decent agreement 
with its counterpart in a laminar viscous disc model. This discrep- 
ancy may arise because in non-isothermal discs, the entropy-related 
corotation torque is associated with singular production of vorten- 
sity on downstream separatrices and we may expect that compared 
with a turbulent disc, such a process is slightly more efficient in a 




Radius 

Figure 9. Upper panel: Entropy profile in turbulent runs with q = 5 X 10~ 6 
and time averaged over 1 500 T or t . For y = 2 X 1 0~ 4 , the entropy profile was 
offset by a factor of 2 for clarity. The dashed line corresponds to the initial 
entropy profile. Lower panel: Same but for the surface density profile. 



Protoplanetary migration in non-isothermal discs with turbulence driven by stochastic forcing 1 1 



laminar viscous disc model. The fair agreement between turbulent 
and laminar torques reveals not only that the total torque experi- 
enced by a protoplanet can be decomposed into a laminar torque 
plus a stochastic torque, which was already suggested by Baruteau 
& Lin (2010) in the case of isothermal disc models, but also that 
the entropy-related horseshoe drag takes similar values in turbu- 
lent and laminar disc models. This suggests that the structure of 
the horseshoe region in turbulent discs should be similar in time 
average to that in laminar disc models, which can be confirmed by 
comparing contour plots of the perturbed entropy for a turbulent 
run to those corresponding to a laminar calculation with similar 
values for the viscous and thermal diffusion coefficients. In Fig.UOl 
we present the perturbed entropy field for a laminar simulation (left 
panel) with k = 4.5 X 10~ 7 and v = 5.4 x 10~ 7 and for a turbulent 
run (middle panel) withy = 6x 10~ 5 (which correspond to effective 
values v ~ 5.4 x 10~ 7 and k ~ 4.5 x 10~ 7 ) at t = 800 T orh . Over- 
plotted are a few streamlines which are indicated with black lines 
and which were computed using the same initial coordinates in the 
R - <p plane for the laminar and the turbulent runs. The right panel 
displays, for the turbulent run, the corresponding entropy field and 
streamlines averaged between 750 and 850 orbits (using field out- 
puts produced every planet orbit) and over 6 different realizations. 
Interestingly, one can see that the averaged streamlines in the tur- 
bulent run have very similar shape to the streamlines in the laminar 
simulation, which is in agreement with the results of Baruteau & 
Lin (2010) in the case of isothermal discs. Here, since the initial en- 
tropy gradient is positive, dynamics in the horseshoe region yields 
a positive entropy perturbation at the inward downstream separatrix 
(0 > (j> p \ R < R p ) whereas a negative entropy perturbation is cre- 
ated at the outward downstream separatrix (0 < (j> p ; R > R p ). It is 
clear that these entropy perturbations for the turbulent run are very 
similar to those corresponding to the laminar calculation, which 
confirms that the thermal diffusion coefficients are indeed close to 
each other in both cases. 



6 DISCUSSION AND CONCLUSION 

In this paper, we have presented the results of 2D hydrodynamical 
simulations which study the interaction of an embedded low-mass 
planet with a turbulent non-isothermal disc. The aim is to investi- 
gate the role of turbulence on Type I migration in non-isothermal 
discs, with particular emphasis on the process of desaturation of the 
entropy-related corotation torque by turbulence. 
We employ the turbulence model of Laughlin et al. (2004) in which 
a turbulent potential corresponding to the superposition of multiple 
wave-like modes is applied to the disc. This turbulence potential 
had been shown previously to reproduce the main statistical proper- 
ties of MHD turbulence in isothermal discs (Baruteau & Lin 2010) . 
In non-isothermal discs, turbulent density and entropy fluctuations 
generate transport of momentum and heat inside the disc and we 
show that both the associated viscous and thermal diffusion coeffi- 
cients scale as y 2 , where y refers to the amplitude of the turbulent 
forcing. This results in a turbulent Prandtl number P r = a R c s H/K 
where a R is the Reynolds stress parameter which is almost con- 
stant with y and we estimate for the first time that P r ~ 0.3 in 
non-isothermal turbulent discs. We caution, however, that such a 
value for the turbulent Prandtl number should be checked against 
3D MHD, non-isothermal simulations in which turbulence is gen- 
erated by the MRI. 

By investigating the evolution of the running time-averaged torque 
experienced by a protoplanet for different values of y, we show that 



turbulence can unsaturate the entropy-related corotation torque. For 
similar values for the viscous and thermal diffusion coefficients, 
we find that the formulae of Paardekooper et al. (2011) for non- 
isothermal Type I migration in laminar disc models reproduce quite 
satisfactorily the time-averaged torques deduced from our turbulent 
runs, provided that the typical amplitude of turbulent fluctuations is 
not too large so that turbulent transport of heat and momentum do 
not significantly modify the background density and temperature 
profiles. For high values of y corresponding to values for the di- 
mensionless a viscosity of Shakura & Sunyaev (1973) a > 10~ 3 , 
we indeed observe discrepancies between the analytical estimate 
and the results of simulations which, to a large extent, are related to 
the alteration of the temperature profile by turbulence. We find also 
that the amplitudes of the running time-averaged turbulent torques 
are slightly smaller than the torques amplitudes deduced from lam- 
inar simulations with initial density and temperature profiles cor- 
responding to the time averaged profiles resulting from turbulent 
runs. This suggests that the entropy-related horseshoe drag in pres- 
ence of turbulence is slightly more saturated compared to an equiv- 
alent laminar disc model with similar vortensity and entropy diffu- 
sion coefficients. 

Although these results were obtained from a disc model in which 
the initial entropy gradient is positive, similar conclusions should 
be drawn in the case of a negative initial entropy gradient. Since 
in that case the amplitude of the total torque is only a small frac- 
tion of the typical amplitude of the stochastic torque, turbulent runs 
would require a much longer convergence time than the simulations 
presented here. For a negative initial entropy gradient, the entropy- 
related corotation torque is positive and can eventually counteract 
the effect of the (negative) differential Lindblad torque, depending 
on the local temperature gradient and provided that the viscous and 
thermal diffusion coefficients are large enough to keep the entropy- 
related corotation torque from saturating. 

The positive results obtained in this paper bring further support 
for the entropy-related corotation torque as a possible solution for 
the discrepancy between the mass-semi-major axis diagram of ob- 
served exoplanets, and that inferred from theories of planetary for- 
mation and evolution (e.g., Ida & Lin 2008, Mordasini et al. 2009, 
Howard et al. 2010, Hellary & Nelson 2012). 
The simulations we have presented have neglected radiative pro- 
cesses and we will present in a future paper simulations with ra- 
diative transfert treated in the flux-limited approximation. We will 
also investigate the effects of turbulence on the dynamics of mul- 
tiple protoplanets located in the vicinity of an opacity transition 
where an entropy jump can act as a planet trap and promotes em- 
bryo growth (Horn et al. 2012). 



ACKNOWLEDGMENTS 

Simulations were performed using HPC resources from GENCI- 
cines (c2012046733). Clement Baruteau is supported by a Herchel- 
Smith Postdoctoral Fellowship of the University of Cambridge. 

REFERENCES 

Balbus, S. A., & Hawley, J. F. 1991, ApJ, 376, 214 
Balbus, S. A., & Papaloizou, J. C. B. 1999, ApJ, 521, 650 
Baruteau, C, & Masset, F. 2008, ApJ, 672, 1054 
Baruteau, C, & Lin, D. N. C. 2010, ApJ, 709, 759 



12 A. Pierens, C. Baruteau and F. Hersant 



Baruteau, C, Fromang, S., Nelson, R. P., & Masset, F. 2011, 

A&A, 533, A84 
Baruteau, C, & Masset, F. 2012. larXiv: 1203.32941 
de Val-Borro, M., Edgar, R. G., Artymowicz, R, et al. 2006, MN- 

RAS, 370, 529 
Fromang, S., & Nelson, R. P. 2006, A&A, 457, 343 
Fromang, S., & Nelson, R. P. 2009, A&A, 496, 597 
Hellary, P., & Nelson, R. P. 2012, MNRAS, 419, 2737 
Howard, A. W., Marcy, G. W., Johnson, J. A., et al. 2010, Science, 

330, 653 

Goldreich, P., & Tremaine, S. 1979, ApJ, 233, 857 
Horn, B., Lyra, W., Mac Low, M.-M., & Sandor, Z. 2012, ApJ, 
750, 34 

Ida, S., & Lin, D. N. C. 2008, ApJ, 673, 487 
Lathrop, D. P., Fineberg, J., & Swinney, H. L. 1992, Phys. Rev. A, 
46, 6390 

Laughlin, G., Steinacker, A., & Adams, F. C. 2004, ApJ, 608, 489 
Li, H, Li, S., Roller, J., et al. 2005, ApJ, 624, 1003 
Lin, D. N. C, & Papaloizou, J. 1986, ApJ, 309, 846 
Masset, F. 2000, A&AS, 141, 165 

Masset, F. S., Morbidelli, A., Crida, A., & Ferreira, J. 2006, ApJ, 
642, 478 

Masset, F. S., & Casoli, J. 2009, ApJ, 703, 857 
Masset, F. S., & Casoli, J. 2010, ApJ, 723, 1393 
Mordasini, C, Alibert, Y., Benz, W., & Naef, D. 2009, A& A, 
501, 1161 

Nelson, R. P., Papaloizou, J. C. B., Masset, F, & Kley, W. 2000, 

MNRAS, 318, 18 
Nelson, R. P. 2005, A&A, 443, 1067 

Ogihara, M., Ida, S., & Morbidelli, A. 2007, Icarus, 188, 522 
Paardekooper, S.-J., & Papaloizou, J. C. B. 2008, A&A, 485, 877 
Paardekooper, S.-J., & Papaloizou, J. C. B. 2009, MNRAS, 394, 
2283 

Paardekooper, S.-J., Baruteau, C, Crida, A., & Kley, W. 2010, 

MNRAS, 401, 1950 
Paardekooper, S.-J., Baruteau, C, & Kley, W. 2011, MNRAS, 

410, 293 

Papaloizou, J. C. B., & Nelson, R. P. 2003, MNRAS, 339, 983 
Pierens, A., Baruteau, C, & Hersant, F. 201 1, A& A, 531, A5 
Ruediger, G., Elstner, D., & Tschaepe, R. 1988, Acta Astron., 38, 
299 

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337 
Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 753 
Tanaka, H, Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257 
Uribe, A. L., Klahr, H, Flock, M., & Henning, T. 201 1, ApJ, 736, 
85 

van Leer, B. 1977, Journal of Computational Physics, 23, 276 
Ward, W. R. 1997, Icarus, 126, 261 



